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Abstract 

Modeling dispersed solid phases in fluids still represents a computational challenge when considering a small-scale coupling in 
wide systems, such as the atmosphere or industrial processes at high Reynolds numbers. A numerical method is here introduced 
for simulating the dynamics of diffusive heavy inertial particles in turbulent flows. The approach is based on the position/velocity 
phase-space particle distribution. The discretization of velocities is inspired from lattice Boltzmann methods and is chosen to match 
discrete displacements between two time steps. For each spatial position, the time evolution of particles momentum is approximated 
by a finite-volume approach. The proposed method is tested for particles experiencing a Stokes viscous drag with a prescribed fluid 
velocity field in one dimension using a random flow, and in two dimensions with the solution to the forced incompressible Navier- 
Stokes equations. Results show good agreements between Lagrangian and Eulerian dynamics for both spatial clustering and the 
dispersion in particle velocities. This demonstrates the suitability of the proposed approach at large Stokes numbers or for situations 
where details of collision processes are important. 
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1. Introduction 

Particle-laden turbulent flows are found in numerous natu¬ 
ral and industrial situations, ranging from droplet growth in 
clouds and dust accretion in early stellar systems, to turbu¬ 
lent mixing in engines and industrial sprays. In such situa¬ 
tions, the processes that need being modeled and quantified in¬ 
volve the fine-scale dynamical properties of the particles, like 
preferential concentration, collisions and coalescences, chemi¬ 
cal reactions, and modulation of the fluid flow by the particles. 
To address specific microphysical issues, one usually study si¬ 
multaneously the turbulent flow and the dispersed phase using 
Eulerian-Lagrangian direct numerical simulation (see, for in¬ 
stance, IBalachandar and Eaton| |2010t |Toschi and Bodenschatz| 
2009 for recent reviews). This approach is particularly suited 
for monitoring the fluctuations occurring at dissipative scales. 
However, direct numerical simulations are quickly too compu¬ 
tationally expensive for studying particle suspensions in real¬ 
istic settings. On the one hand, a large-scale system, such as 
a chemical reactor, an atmospheric cloud or a protoplanetary 
disk, contains a prohibitively large number of particles. On the 
other hand, the substantial Reynolds numbers of natural and 
industrial flows require the use of large-scale models, such as 
large-eddy simulations. Eulerian-Lagrangian methods, where 
the dispersed phase is modeled by point particles, show some 
advantages: they allow for an easy implementation of polydis- 
persity and are rather insensitive to subgrid-scale fluctuations. 
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at least for particles with a large-enough response time (Wang 


and Squire^ 1 1996| ). For particles with smaller inertia, one relies 


on the use of stochastic Langevin models ( [Pozorski and Apte] 
2QQ9[ [Shotorban and Mashayek| |2006| ). In addition the con¬ 


straints on the number of particles can be relaxed using super¬ 
particles approaches, which then necessitates some modeling 
for collisions ( |Rein et ar||2010|[Shima et ari|2009| ). However, 
as stressed for instance by Portela and Oliemans ( |2006| ), La¬ 
grangian methods prove some difficulties in correctly predict¬ 
ing modifications of the carrier flow by the dispersed phase, 
particle-to-particle interactions and near-wall effects. 

Some of these shortcomings can be circumvent using 
Eulerian-Eulerian methods (see |Fox[|2012[ for a review). The 
main difficulty then relies in finding a fluid description of the 
dispersed particulate phase. In principle, this is achieved by 
prescribing a closure for the kinetic hierarchy of moment equa¬ 
tions. When considering an ensemble of particle trajectories 
(jCp, Vp), one naturally introduces the phase-space density 


fix, V, t) = {six^it) - x)^(Vp(0 - v)), (1) 

where the fluid velocity realization is fixed and the average (•) 
is both over the particle ensemble (different realizations of the 
initial conditions and/or average over a large number of par¬ 
ticles) and over the molecular diffusion of the particles (with 
diffusion constant k). The phase-space density then solves the 
diffusion-transport equation 


d,f + v- V,/ - V, • [/ Tix, V, 0/mp] - V, • (^V,/) = 0, (2) 

where T is the force exerted by the fluid on a particle located 
at X with a velocity v and is the particle mass. This Fokker- 
Planck equation is exact and fully describes the dynamics of 
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small particles in the phase space. The drift terms are com¬ 
pletely prescribed by a given realization of the fluid flow. To ob¬ 
tain an Eulerian description of particles dynamics that depends 
on the spatial variable x only, the usual approach consists in 
deriving the equations for the various moments of the velocity 
V. To close the resulting hierarchy, additional assumptions are 
needed. They naturally arise when focusing on given asymp¬ 
totics (see, e.g., |Carrillo et al. 2008| ). For instance, when the 
particles experience a very strong viscous drag with the flow 
(small Stokes numbers), an effective particle velocity can be 
written { Maxey | 1987| ) leading to close this hierarchy at the 
zeroth order and to write an explicit equation for the trans¬ 
port of particle density. This then serves as a basis for deriv¬ 
ing subgrid-scale models for large-eddy simulations (see, e.g., 
[Shotorban and Balachandar||2007| ). First-order closures lead to 
writing an evolution equation for a particle velocity held that is 
coupled to the fluid flow. Again, such methods are limited to 
asymptotically small values of the Stokes number, as they are 
inadequate to deal with multi-streamed particle distributions. 
It is indeed well known that the trajectories of finite-Stokes- 
number particles can cross, leading to the formation of regions 
where the particle velocities are multivalued and cannot be de¬ 
scribed in terms of a spatial held. This phenomenon is usually 
referred to as particle-trajectory crossing ( |Chen[ [20001 ), sling 
ejfect ( [Falkovich et al.[ |2002| ) or caustic formation ( [Wilkinson| 
land Mehhg 2005 ) and has important impacts in estimating col¬ 
lision rates (see, e.g.. Bee et al.| 2005). Higher-order closures, 
such as ten-moment equations, can also be proposed depending 
on the specific forces applied on the particles. They account 
for the dispersion in particle velocities and can thus catch some 
aspects of multi-streaming (using either algebraic or quadrature 
closures Desjardins et al. 2008[ Fevrier et~5H|2005| ). 

In principle, accessing the full multi-streaming dynamics of 
particles requires solving the kinetic equation 0 in the entire 
(2 X J)-dimensional position-velocity phase space. (J denotes 
here the dimension of the physical space.) A clear difficulty 
which is then faced is the prohibitive computational cost of in¬ 
tegrating a partial differential equation in a space with such a 
large dimensionality. Attempts have nevertheless been made by 
decreasing the number of relevant degrees of freedom. This can 
be easily done, for instance, by considering one-dimensional 
flows (see, e.g., Carrillo et al. 2008| and references therein). 
Other approaches are based on the physical observation that the 
distribution of particle velocities is usually rather concentrated 
along a given number of branches in phase space. This led for 
instance p!^iu et al.| ( 2Qll| ) to capture implicitly the velocity dis¬ 
persion by applying a level-set method in phase space. The 
efficient implementation of this procedure to high-dimensional 
turbulent situations still represents a real challenge. 

Here, we propose an alternative approach that consists in de¬ 
grading the resolution in velocities and to apply computation¬ 
ally efficient ideas inspired from Lattice-Boltzmann methods 
( |Chen and Doolen[|l998| ). The discrete values of velocities are 
chosen to correspond exactly to discrete displacements between 
two time steps on a fixed spatial lattice. The time evolution of 
fix, v,t) is then approximated by splitting the spatial advection 
on the lattice and the acceleration of particles, which is inte¬ 


grated using a finite-volume scheme. This gives access to a 
full phase-space particle distribution that naturally catch multi¬ 
streaming. This method is relevant to cases where diffusion is 
responsible for a broadening of the particle velocity dispersion 
and it applies to any kind of force ‘F acting on the particles. Af¬ 
ter describing the algorithm in Sec.|^ we present some qualita¬ 
tive and quantitative tests for very small heavy particles whose 
dynamics is dominated by diffusion and a Stokes drag with the 
fluid velocity. The proposed lattice-particle method is directly 
compared to the results of direct Lagrangian simulations. Sec¬ 
tion |3] is devoted to the one-dimensional case with a random 
Gaussian flow. Section]^ shows the results of coupled Navier- 
Stokes and lattice particles simulations in two dimensions for 
turbulent flows, either in the direct cascade of enstrophy or in 
the inverse cascade of kinetic energy. 

2. Description of the method 

The solutions fix, v,t) to the Liouville (or Fokker-Planck) 
equation ^ are defined in the full position-velocity phase- 
space Q X where Q designates a J-dimensional bounded 
spatial domain. To simulate numerically the dynamics, we di¬ 
vide the phase-space in (2 x J)-dimensional hypercubes. We 
then approximate fix, v,t) as a piecewise-constant scalar held 
on this lattice. Positions are discretized on a uniform grid with 
spacing Ax in all directions. In principle, / has to be defined 
for arbitrary large velocities. We however assume that relevant 
values of v are restricted to a bounded interval [-Vmax, VmaxY . 
where V^ax has to be specified from physical arguments based 
on the forces F applied on the particles. Velocities are assumed 
to take values, so that the grid spacing reads Av = 2 V^axINv- 
Figure illustrates the phase-space discretization in the one¬ 
dimensional case with Ny = 5. The various cells in position- 
velocity contain a given mass of particles. All these particles 
are assumed to have a position and velocity equal to that at the 
center of the cell. 



Figure 1: Sketch of the lattice dynamics in the (x, v) position-velocity phase 
space. The solid curve is the fluid velocity profile; the greyscale tiling repre¬ 
sents the discretization of particles mass in phase space. The black horizontal 
arrows show advection, while the blue and red vertical arrows are forcing and 
diffusion, respectively. 


The three phase-space differential operators appearing in 
Eq. g., namely the advection, the particle forcing, and the dif¬ 
fusion, are applied one after the other, following an operator 
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splitting method ( [LeVeque] 2002| ). For the advection step, we 
use a technique inspired from the Lattice-Boltzmann method 
(see, e.g., |Succi[[20QT] ). The time stepping is chosen, so that a 
discrete velocity exactly matches a shift in positions by an in¬ 
teger number of gridpoints. Namely, we prescribe Ax = AvAt. 
All the particle phase-space mass located in [-Av/2, Av/2] does 
not move; that in [Av/2, 3Av/2] is shifted by one spatial grid- 
point to the right and that in [-3Av/2, -Av/2] to the left, etc. 
All the mass is displaced from one cell to another according 
to its own discrete velocity value. This evolution is sketched 
by black horizontal arrows in Fig[2 This specific choice for the 
time-stepping implies that the advection (in space) is treated ex¬ 
actly for the discrete system. The next steps consist in applying 
the force acting on the particles and the diffusion. The corre¬ 
sponding terms in Eq. ^ are conservation laws, which suggests 
using a finite-volume approximation. The time evolutions due 
to forcing and diffusion are performed successively. In both 
cases, we use classical schemes (see below), where zero-flux 
conditions are imposed on the boundary of [-Vmax, Fmax]^- The 
force is evaluated using the values of v at the centers of the cells 
and Vv/ is approximated using finite differences. These steps 
are illustrated by the horizontal blue and red arrows in Fig.[2 

A few comments on the convergence and stability of the pro¬ 
posed method. Clearly, except for specific singular forcings, 
all the linear differential operators involved in ^ are expected 
to be bounded^ We can thus invoke the equivalence (or Lax- 
Richtmyer) theorem for linear differential equations that en¬ 
sures convergence, provided the scheme is stable and consistent 
(see |LeVeque[ |2002| ). 

For the operator associated to particle acceleration, we use in 
this study either a first-order upwind finite-volume scheme or 
a higher-order flux limiter by following the strategy proposed 
by |Hundsdorfer et ^ ( 1995| ). The upwind scheme is first-order 
accurate and is well-known for being consistent and stable if it 
satisfies the Courant-Friedrichs-Lewy (CFL) condition. This 
requires that the time needed to accelerate particles by the grid 
size Av is larger than the time step At. This leads to the stability 
condition 


CFL = {AtjAv) max |‘F(jc, V, Ol/^p < 1- (3) 

x,v,t 

The upwind scheme is however known to suffer from numeri¬ 
cal diffusion, and obviously, one should only expect to recover 
the correct dynamics only when the numerical diffusion /Cnum is 
much smaller than the physical one k. The average numerical 
diffusion can be estimated as (/Cnum) ~ {T)AtjAv. To limit the 
effects of this numerical diffusion, we have also used a flux- 
limiter scheme. While taking benefit of a higher-order approxi¬ 
mation where the field is smooth, it uses the ratio between con¬ 
secutive flux gradients to reduce the order in the presence of 
strong gradients only. The limiter is a nonlinear function of the 
phase-space density field and the stability is ensured provided 
that it is total-variation diminishing (TVD; see |LeVeque||2002| ). 


^Notice that, although the velocity might explicitly appear in the force T, 
we only solve for a compact domain of velocities, thus preventing divergences. 


Among the various TVD limiters available on the market, we 
choose the scheme of |Koren| ( |1993| ) with parameter 2/3. 

For the term associated to diffusion, the flux at the interface 
between two velocity cells is computed using finite differences. 
The resulting finite-volume scheme is thus equivalent to com¬ 
pute a discrete Laplacian on the velocity mesh. The stability 
condition is then given by 


KAt 1 
Av^ 2 


(4) 


To summarize, the stability and convergence of the proposed 
method is ensured when both ^ and (|^ are satisfied. 

From now on we restrict ourselves to small and heavy par¬ 
ticles whose interaction with the carrier fluid is dominated by 
viscous drag and diffusion. In that case, we have 


dvp 

dt 



(vp -M(A:p,0) + ^f^T](t), 


(5) 


where T](t) is the standard J-dimensional white noise and the 
fluid velocity field u(x, t) is prescribed and assumed to be in a 
(statistically) stationary state. This Stokes drag involves the vis¬ 
cous particle response time Tp = 2p^a^ 1(9pfv), with a the diam¬ 
eter of the particles, v the viscosity of the fluid, Pp » pf the par¬ 
ticle and fluid mass densities, respectively. Inertia is quantified 
by the Stokes number St = Tp/rf, where Tf is a characteristic 
time of the carrier flow. The diffusion results from the random 
collisions between the considered macroscopic particle and the 
molecules of the underlying gas. Assuming thermodynamic 
equilibrium, the diffusion coefficient reads k = Ik^Tl(m^ Tp), 
where kp is the Boltzmann constant and T the absolute tem¬ 
perature. The effect of diffusion is measured by the non- 
dimensional number = KTijUj (Uf being a characteristic 
velocity of the fluid flow). 

Such a specific dynamics leads to appropriate estimates for 
the bound V^ax in particle velocity. One can indeed easily check 
that when /c = 0, we always have |Vp| < max^ Jm(jc, 01- In a 
deterministic fluid flow, as for instance when u is stationary, 
this gives the natural choice Vmax = maxjj: ^ |m(jc, 01- However, 
in most situations, the maximal fluid velocity is not known a 
priori. One then relies on the statistical properties of m, as for 
instance its root-mean square value Usually 

the one-time, one-point statistics of fluctuating velocity fields 
(being random or turbulent) are well described by a Gaussian 
distribution. This ensures that by choosing Vj^ax = 3 Wrms. the 
probability that a particle has a velocity out of the prescribed 
bounds is less than 1%. Such estimates are rather rough. In 
practice, it is known that the typical particle velocity decreases 
as a function of the Stokes number. It was for instance shown by 
) that (|Vpp) oc u^^jSt at very large Stokes 
numbers. An efficient choice for Vj^ax should account for that. 

In the next two sections we investigate two different cases: 
First a one-dimensional random Gaussian carrier flow with a 
prescribed correlation time and, second, a two-dimensional tur¬ 
bulent carrier flow that is a solution to the forced incompressible 
Navier-Stokes equation. 
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3. Application to a one-dimensional random flow 


3.1. Particle dynamics for d = \ 

In this section, we test our method in a one-dimensional sit¬ 
uation. For that, we assume that the fluid velocity is a Gaussian 
random fleld, which consists in the superposition of two modes 
whose amplitudes are Ornstein-Uhlenbeck processes, namely 


u{x, t) 
dA,(0 
dt 


A\{f) cos(27rx/L) -l- A 2 {f) sin(2;7rx/L) 


- Aft) + 

Tf 



( 6 ) 

(7) 


where the are independent white noises with correlations 
{fi{t)fi{t')) = u^^6{t - t'). This flow is by deflnition fully 
compressible (potential) and spatially periodic with period L. 
It is characterized by its amplitude {{u{x, = Wnns and its 

correlation time Tf, which are flxed parameters. We focus on 
the case when the Kubo number Ku = Tf Wrms/7^ is of the order 
of unity. 

We next consider particles suspended in this flow and fol¬ 
lowing the dynamics The relevant Stokes number is then 
deflned as St = Tp u^msIL and the relative impact of diffusion 
is measured by = kLIu^^. When diffusion is neglected 
(7C ^ 0), the particles distribute on a dynamical attractor (see 
Fig|^ Left) whose properties depend strongly on St. These 




X X 

Figure 2: Instantaneous snapshots of the particle positions in the (x, v) plane for 
St ^ 2 and with no diffusion (TC = 0, Left) and with diffusion CTC = 6 x 10“^, 
Right). The folded structures are spread out by diffusion. 




Strange attractors are typically fractal objects in the phase space 
and they are characterized by their fractal dimension spectrum 
( Hentschel and Procaccia[p^8^ . The locations of particles are 
obtained by projecting theses sets on the position space and 
might thus inherit the associated clustering ( |Bec 2003| ). The 
dimension that is relevant for binary interactions between parti¬ 
cles is the correlation dimension D 2 , which relates to the prob¬ 
ability of having two particles within a given distance, namely 

P 2 (r) = IP(l4’- -^^(01 <r)~ {rlLf\ for r L, (8) 


where and x® denote the positions of two different parti- 
cles. Note that as we consider u to be in a statistically stationary 
state, P 2 is independent of time . The correlation dimension D 2 
varies from £>2 = 0 for a point concentrations to D 2 = 1 for a 
homogenous mass distribution . In the example of Fig. (Left) 
D 2 ~ 0.7. The variations of D 2 as a function of the Stokes 
number are displayed in the inset of Fig. D 2 indeed varies 
from 0 at small Stokes numbers to values close to one. For 
St = 0, the particles concentrate on a point; their distribution 
is said to be atomic and D 2 = 0. This is a consequence of the 


compressibility of the one-dimensional (potential) flow. Ac¬ 
tually this behavior persists for flnite Stokes numbers, up to a 
critical value St^, as shown by Wilkinson and Mehlig (2003) in 
the case where Tf L/wnns (that is Ku ^ {)). We observe here 
St^ ^ 0.6. For St > St*, the dimension increases and tends to a 
homogenous distribution (D 2 = 1) at large particle inertia. 

When one has only access to the Eulerian density of particles, 
the distribution of distances cannot be directly inferred from 
(|8J. One then relies on the coarse-grained density of particles 


pr(x, 0 = J" hv' J" dv f(x -h x\v,t). 


(9) 


It is known that, under some assumptions on the ergodicity of 



the same as the radial distribution function. This quantity will 
be used in the next sections to address the physical relevance of 
the lattice-particle method. It is of particular interest when con¬ 
sidering collisions between particles. Indeed, as explained for 
instance by [Sundaram and Collins | \\991) , the ghost-collision 
approximation leads to write the collision rate between particles 
as the product of two contributions: one coming from cluster¬ 
ing and entailed in the radial distribution function, and another 
related to the typical velocity differences between particles at 
a given distance. This second quantity relates to the particle 
velocity (first-order) structure function 


5i(r) = (|v[,' 


,( 1 ) 




•vn||x^ 


( 1 ) 


..( 2 ) 


\-r)^ 


( 10 ) 


This is the average of the amplitude of the velocity difference 
between two particles that are at a given distance r. As the 
probability of distances, this q uantity behaves as a power law 


Si(r) ~ for r L (see, e.g.,|Bec et al. 2005). The exponent 


fi, shown in the inset of Fig. [^decreases from 1 Sit St = 0, cor¬ 
responding to a differentiable particle velocity fleld, to 0 when 
St ^ 00 , which indicates that particle velocity differences be¬ 
come uncorrelated with their distances. Again, when working 
with the phase-space density one cannot use ( p^ but relies on 


^i(r) = 


{Jdvfdv'f(x,v)f(x + ry)\v-v'\) 


( 11 ) 


As the second-order moment of the coarse-grained density, this 
quantity will also be used as a physical observable for bench¬ 
marking the method. 

In the above discussion, we have neglected the effects of dif¬ 
fusion. It is for instance expected to alter clustering properties 
by blurring the particle distribution at small scales. This is il¬ 
lustrated in Fig. where one can compare the instantaneous 
phase-space particle positions in the absence of diffusion (Left) 
and when it is present (Right) at the same time and for the same 
realization of the fluid velocity. At large scales, identical pat¬ 
terns are present, but diffusion acts at small scale and smoothes 
out the fine fractal structure of the distribution. One can easily 
estimate the scales at which this crossover occurs. Diffusion is 
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responsible for a dispersion in velocities that can be obtained 
by balancing Stokes drag and diffusion in the particle dynam¬ 
ics, namely v^/Tp k , so that Vd ~ This dispersion 

in velocity is responsible for a dispersion in positions on scales 
of the order of Id = Tp Vd ~ Hence, 

when diffusion is small enough and ^d ^ L, the spatial distri¬ 
bution of particles is unchanged by diffusion at length scales 
r » ^d. and the probability that two particles are at a distance 
less than r behaves as p^ir) ~ (rlL)^^. For r diffusion 

becomes dominant, the particles distribute in a homogeneous 
manner and oc with d = I being the space dimension. 
By continuity at r = ^d^ we get p^ir) ~ at small 

scales. 



Figure 3: Cumulative probability pj(r) of inter-particle distances for various 
diffusivities k and for St ^ 2. One observes at low diffusivities and for r > 
a behaviour oc (rlL)^^ with D 2 ~ 0.1 < d = 1, followed at small scales by 
uniform particle distribution where pj(r) oc r. As k increases, the transition 
is moved to larger values of r. The vertical dashed represent in each case the 
estimate £d ~ for this transition. Inset: variations of the correlation 

dimension D 2 and of the scaling exponent of the first order particle structure 
function, as a function of the Stokes number St in the case of the random fluid 
velocity defined by Eqs. 

This picture is confirmed numerically as shown in Fig. 
which represents the scale-behavior of p^ir) for a fixed Stokes 
number and various values of the diffusivity k. One clearly ob¬ 
serves the homogenous distribution cx at small scales and the 
fractal scaling cx intermediate range. The predicted 

transition between the two behaviors is indicated by the verti¬ 
cal lines at the diffusive scale ^d- A homogeneous distribution 
is recovered for r < ^d/10- 

Velocity statistics are also altered by the presence of diffu¬ 
sion. The structure function Si(r) is expected to behave as 
for r L and to saturate to a constant value when r ^d- 
By continuity, the value of this plateau should be ~ ~ 

Note finally that the slow convergence cx as ^ 0 
implies that very small values of the diffusion are needed in or¬ 
der to clearly recover the statistics of diffusive-less particles as 
an intermediate asymptotics. 

3.2. Lattice-particle simulations 

We now turn to the application of the lattice-particle method 
described in Sec. |2] to this one-dimensional situation. We 


compare the results to Lagrangian simulations where we track 
the time evolution of Np particles randomly seeded in space 
with zero initial velocity. We choose and normalize the ini¬ 
tial phase-space density /(x, v, 0) to match the Lagrangian set¬ 
tings. The distribution is uniform over the cells, concen¬ 
trated on a vanishing velocity and the total mass is such that 
Yjij f(Xi,Vj,t) Ax Av = Np. In all simulations, the maximum 
velocity is set to Fmax = ^rms = 1 and we have chosen L = In 
and Tf = 1. In these units, the time step is kept fixed at 
At = Ax/Av = 2“^7r 0.05. The number of discrete veloci¬ 

ties is of the form 2" -r 1 and is varied between = 3 to 129. 
The number of spatial collocation points is then given by 2”^^ 
and thus varies between Nx = 128 to 8192. Note that, because 
of the CFL condition ([^, this choice restricts the number of 
discrete velocities that can be used to Vv < 1 + 128 St. 



0 TT ItT 

X 

Figure 4: Position-velocity phase-space positions of Lagrangian particles 
(black dots) on the top of the field obtained by the lattice-particle method (col¬ 
ored background). The diffusivity is here k = 10“^ (7C ^ 6.28.10“^) and St ^ 2. 

Figure [^represents simultaneously the phase-space distribu¬ 
tion of Lagrangian particles and the numerical approximation 
obtained by the particle-lattice method for Ny = 129. Clearly, 
one observes that the model fairly reproduces the distribution 
of particles, including the depleted zones, as well as the more 
concentrated regions. Furthermore, the method is able to catch 
multivalued particle velocities. We have for instance up to three 
branches in v for x ^ 37r/2. It is important to emphasize that 
numerical diffusion is of course present, and that it has to be 
smaller than the physical diffusion k in order for the model to 
be consistent with the Lagrangian dynamics. 

To get a more quantitative insight on the convergence of the 
method, we next compare the coarse-grained densities obtained 
from the Lagrangian simulation and the lattice-particle approx¬ 
imation of the phase-space density. The first, denoted pjr is 
computed by counting the number of particles contained in the 
different boxes of a tiling of size r. The second is written as 
and is obtained by summing over velocities and coarse-graining 
over a scale r the phase-space density obtained numerically. To 
confirm the convergence of the method, we measure for a fixed 
r the behavior of the X 2 -norm of the difference between p|r and 
pf, namely 

llPr -PrII = , (12) 
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where the angular brackets (•) encompass a spatial and a time 
average. Figureshows the behavior of the relative X 2 -error as 



Ny 

Figure 5: Relative X 2 -error of the lattice-particle method for evaluating the 
coarse-grained density pr over a scale r = L/128 as a function of the number of 
velocity gridpoints Ny and for various values of the Stokes number, as labeled. 


a function of the number of velocity gridpoints for various 
values of the Stokes number St and for a given scale r. One 
observes that the error decreases when the resolution increases, 
giving strong evidence of the convergence of the method. The 
error is found proportional to the velocity grid spacing Av, indi¬ 
cating that the method is first order. The constant is a decreasing 
function of the Stokes number. This indicates that the method 
is more accurate for particles with strong inertia. The reason 
for this trend will be addressed in the sequel. 

To assess the ability of the proposed method to reproduce 
physically relevant quantities, we now compare statistics ob¬ 
tained using the lattice method with those using a Lagrangian 
approach. We focus on the clustering and velocity difference 


properties that were introduced and discussed in Section 3.1 



r/L 

Figure 6 : Convergence of the second-order moment of the coarse-grained 
density <p^>, which is shown as a function of r for St ^ 1.9, TC = tt 10“^, and 
various lattice velocity resolutions Ny, as labeled. 

Figure shows for given values of the Stokes number and of 
the diffusivity, the second-order moment of the coarse-grained 
density as a function of r and various values of the res¬ 

olution in velocities, together with the value ((p|r)^) obtained 
with 10^ Lagrangian particles. One observes that the curves 
approach the limiting behavior from below when the number 


of gridpoints becomes larger (i.e. when Av ^ 0). At suffi¬ 
ciently high velocity resolutions, the method is able to capture 
the large-scale properties of the concentration of the particles. 
The second-order moment of density then saturates to a value 
lower than that expected from Lagrangian measurements. The 
situation is very different at very low resolutions where the data 
obtained from the lattice-particle method deviates much, even 
at large scales. This corresponds to the case when the numerical 
diffusion in velocity is larger than the physical diffusion. 

These strong deviations stem from a non-trivial effect of dif¬ 
fusion that lead to finite-scale divergences of the solutions asso¬ 
ciated to different values of In the absence of diffusion, there 
is a finite probability that an order-one fraction of mass gets 
concentrated on an arbitrary small subdomain of the position- 
velocity phase space. This corresponds to a violent fluctuation 
where the local dimension approaches zero. At the time when 
this occurs, the mass distribution associated to a finite value 
of the diffusion will get stacked at a scale Because of the 
chaotic nature of the particle dynamics, the two mass distribu¬ 
tions, with and without diffusion will experience very differ¬ 
ent evolutions and diverge exponentially fast. Such a strong 



Figure 7: Three snapshots of the Lagrangian particles (black dots, for = 0) 
and of the lattice-particle Eulerian solution (colored background) in the (x, v) 
plane for different times: At r < 4 (Left) the solution is well approximated at 
large scales; At r = r.*. (Center) an order-unity fraction of the mass is concen¬ 
trated on a scale less than AX t > (Right), the Eulerian and Lagrangian 
solutions diverge exponentially fast with differences appearing at the largest 
scales.. 

clustering event followed by the divergence of the solutions, 
is shown in Fig. [7] Starting from a correctly reproduced dis¬ 
tribution, the major part of non-diffusive Lagrangian particles 
concentrate into a subgrid region while the Eulerian approxi¬ 
mation is stacked at scales of the order of At a later time, 
the two distributions diverge and the diffusive particles fill faster 
larger scales. The probability with which one encounters such 
a configuration strongly depends on the Stokes number and on 
the spatial dimension. In the one-dimensional case, such events 
are rather frequent but become sparser when the Stokes number 
increases. This is essentially due to the compressibility of the 
carrier flow. For incompressible fluids in higher dimensions, 
we expect a negligible contribution from these events. 

To close this section on one-dimensional benchmarks of the 
lattice-particle method, we report some results on velocity dif¬ 
ference statistics. For that, we have measured the first-order 
structure function Si{r) of the particle velocity, using ( p^ in the 
Lagrangian case and ([TTlfor solutions obtained with the lattice- 
particle method. Figure]^ shows the relative error of Si{r) for 
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fixed values of the separation r, the Stokes number, and the dif- 
fusivity, as a function of the velocity resolution. Clearly, when 
the number of gridpoints Ny increases, the error decreases, fol¬ 
lowing a law approximatively proportional to the grid spacing 
Av. The inset shows the same quantity but, this time, for a 
fixed resolution (Ny = 33) and as a function of the Stokes num¬ 
ber. One clearly observes a trend for this error to decrease with 
St. There are two explanations for this behavior. First, as seen 
above, there are strong clustering events leading to differences 
between the Lagrangian and lattice solutions that can persist for 
a finite time. When the Stokes number increases, such events 
become less probable. The second explanation relies on the fact 
that particles with a larger Stokes number experience weaker 
velocity fluctuations. This implies that for a fixed value of Vmax, 
the particle velocity is more likely to be fully resolved at large 
values of St. As seen in the inset of Fig. the downtrend of 
the error is compatible with a behavior oc Sr^^^. It might thus 
be proportional to the expected value of the root-mean-squared 
particle velocity when » 1 ( [Abrahamson 1975), favoring 
the second explanation. We will turn back in Sec. |4.2.1| for the 
two-dimensional case on the effect of the maximal resolved ve¬ 
locity Vmax onto the convergence of the lattice-particle method. 



Figure 8: Relative error between the particle velocity structure function S^(r) 
obtained from the lattice-particle method and that from Lagrangian aver¬ 
ages, as a function of the number Ny of velocity gridpoints. Here, the Stokes 
number is fixed Xr ^ 1.9 and TC = 2n 10 “^. Inset: same quantity but for//v = 33 
and as a function of the Stokes number St. 


4. Application to incompressible two-dimensional flows 

We extend in this section the lattice-particle method to two- 
dimensional flows. For the particle acceleration, we again make 
use of a flux-limiter scheme. 

4.1. Cellular flow 

We first consider a fluid flow that is a stationary solution to 
the incompressible Euler equations (and to the forced Navier- 
Stokes equations). It consists of a cellular flow field, a model 
that have often been used to investigate mixing properties, 
as well as the settling of heavy inertial particles (see, e.g., 
|Bergougnoux et al.]|2Q14t[Maxey and Corrsin|[l986| ). The ve¬ 
locity field is the orthogonal gradient of the L-periodic bimodal 
stream function y) = U sin(7r(v -f y)/L) sin(7r(x - y)/L) (the 



Figure 9: Particles stationary distribution inside a tilted cellular flow along with 
the density field from the lattice method. The value of the diffusivity is TC = 
Stt X 10“^. Left: St = l/(2n). Right: St = I/tt. These simulations were 
performed on a lattice with 1024^ spatial gridpoints associated to 19^ discrete 
velocities. The distributions have been here spatially shifted in order to avoid 
having the concentration point (0,0) at the origin. 


typical velocity strength is here denoted by U). The cellular 
flow has been here tilted by an angle njA in order to avoid any 
alignment of the separatrices between cells with the lattice that 
leads to spurious anisotropic effects. 

Figure shows two snapshots for two different values of 
St = T^UIL of the stationary particle distribution (black dots), 
together with the density field evolved by the lattice-particle 
dynamics. For the smallest Stokes number (Left panel), one 
observes that the particle distribution is concentrated along the 
separatrices between the different cells. One also observes that 
it develops entangled structures in the vicinity of the hyperbolic 
stagnation points of the flow. These loops, which are aligned 
with the stable direction, corresponds to oscillations in the par¬ 
ticle dynamics that occurs when their inertia makes them cross 
the unstable manifold with a too large velocity. At larger St, 
the particle distribution is somewhat broader but is this time 
centered on specific trajectories that do not perform the afore¬ 
mentioned oscillation but rather cross ballistically the hetero¬ 
clinic separatrices. These results show that the lattice-particle 
method is able to reproduce the complex dynamics of particles 
in a two-dimensional steady flow. The fine structures of the spa¬ 
tial distribution are fairly reproduced, as long as the numerical 
diffusion surpasses numerical errors. 

4.2. Heavy particles in 2D turbulence 

We next turn to the study of the model in non-stationary fluid 
flows that are solutions to the forced two-dimensional incom¬ 
pressible Navier-Stokes equations 

dfU -h M • Vm = -S/p -h vV^M - au + f, V • m = 0. (13) 

The linear damping term involving the coefficient a originates 
from Ekman friction (in geophysical flows), Rayleigh friction 
(in stratified fluids) or the friction induced by the surrounding 
air in soap-film experiment. The flow is maintained in a sta¬ 
tistical steady state by the forcing / that is assumed to be con¬ 
centrated over a specific scale if. The fluid velocity field u is 
computed numerically using a pseudo-spectral, fully de-aliased 
GPU solver for the vorticity-streamfunction formulation of the 
Navier-Stokes equation ( p3] ). 


7 

















The two-dimensional Navier-Stokes equation is known to 
develop two cascades (see, e.g., Botfetta and Ecke| 20121. Ki¬ 
netic energy undergoes an inverse cascade with a rate s toward 
the large scales ^ where it is dissipated by the linear damp¬ 


ing. The enstrophy {clt) (where co 


Jx^y 


Jyl/lX 


designates 


the vorticity) experiences a direct cascade to the small scales 
I if with a rate rj and is then dissipated by molecular vis¬ 
cosity. These different cascades are associated to different be¬ 
haviors of the velocity power spectrum. For k «: kf oc 
the inverse energy cascade promotes a Kolmogorov law, 
as in the three-dimensional direct cascade. At small scales, i.e. 
for k » kf in the direct enstrophy cascade, the flow is charac¬ 
terized by long-living vortices and the spectrum follows a k~^ 
Batchelor-Kraichnan law with a logarithmic correction. 

Dimensional analysis predicts that the direct enstrophy cas¬ 
cade is associated to a unique timescale tq = Inves¬ 

tigating heavy particle dynamics at the small scales of two- 
dimensional turbulence thus requires comparing their response 
time to tq. The relevant parameter is then the Stokes number 
defined as St = t^Itq. For St 1, particles almost follow 
the flow and tend to distribute homogeneously in space. When 
5^ » 1, they completely detach from the fluid and experience 
a ballistic motion leading again to a space-filling distribution. 
Non-trivial clustering effects occur when the Stokes number is 
order one. This is illustrated in Fig. [TO which shows a snapshot 



Figure 10: Snapshot of the position of particles (black dots) for St = 0.1. The 
colored background shows the vorticity field obtained from a 1024^ direct nu¬ 
merical simulation with a large-scale forcing at wavenumbers 1 < |A:| < 4. 

of the particle distribution in the position space on top of the 
turbulent vorticity field in the direct enstrophy cascade. Due to 
their inertia, particles are ejected from vortices and concentrate 
in high-strain regions. There, the combination of stretching, 
folding and dissipation induced by their dynamics makes them 
converge to a dynamical attractor with fractal properties. Such 
a behavior is quantitatively measured by the correlation dimen¬ 
sion D 2 defined in Eq. The evaluation of D 2 as a func¬ 


tion of St resulting from Lagrangian simulations is presented in 
Fig. 0 At St = 0, unlike in the one-dimensional case where 
the dimension of the attractor is 0, particles follow the stream¬ 
lines of the incompressible two-dimensional flow, fill the posi¬ 
tion space, and hence D 2 = 2. Clustering then increases with 
inertia to attain a minimum at St ^ 0.2. It then decreases again 
as the velocity of particles separate from that of the fluid and 
disperse in the velocity space, leading to a space-filling distri¬ 
bution D 2 = 2 when St ^ 00 . 



Figure 11: Correlation dimension D 2 and scaling exponent of the particle 
velocity first-order structure function as a function of St in the two-dimensional 
direct enstrophy cascade. 


The velocity distribution of particles is itself having a behav¬ 
ior that is very similar to the one-dimensional case. This is clear 
from Fig. [TT] where the scaling exponent of its first-order 


structure function (see Eq. ( p^ ) is represented as a function of 
the Stokes number. For S'f 1, the particles are as if advected 
by a smooth velocity field and ~ 1. When St > 1, particles 
with very different velocities can come arbitrarily close to each 
other and ^ 0. 

Particle properties in the inverse energy cascade are more 
difficult to characterize because of the scale-invariance of the 
velocity field. In particular, neither the moments of the coarse¬ 
grained density nor the particle velocity structure functions dis¬ 
play any scaling behavior. What has been nevertheless observed 
numerically by |Boffetta et al.| ( |2004| ) is that the particle spatial 
distribution is dominated by the presence of voids whose sizes 
obey a universal scaling law. Chen et al. ( 2006| ) argued that 
such voids are related to the excited regions of the flow and that 
particles tend to follow the calm regions where the zeros of the 
fluid acceleration are more probable. 

In the sequel we apply the lattice method to both the di¬ 
rect and the inverse two-dimensional cascades. Resolving both 
cascades in the same simulation would require a tremendous 


scale separation and thus number of gridpoints (see Boffetta and 


[Musacchiol |2010| ). For that reason we consider the two cases 
separately. 


4.2.1. Direct enstrophy cascade 

The fluid flow is integrated by a pseudo-spectral method on 
a uniform square spatial grid using a stream-function formu- 
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lation of the Navier-Stokes equation ( p^ . To maintain a de¬ 
veloped turbulent state, a stochastic forcing is applied in the 
wavenumber shell 1 < \k\ <4 of Fourier space while the kinetic 
energy accumulating at large scales is removed by a linear fric¬ 
tion. The particle dynamics is simulated using a spatial lattice 
with the same resolution as the fluid and with various numbers 
of discrete velocities. The acceleration step is done via a 
flux limiter scheme as described in Sec. [21 Results are com¬ 
pared to particle trajectories obtained from Lagrangian simula¬ 
tions. Figure \V7\ shows the instantaneous particle distributions 
obtained from the two approaches. The qualitative agreement is 
excellent, reproducing correctly depleted zone as well as more 
concentrated regions. 



Figure 12: Snapshot of the position of Lagrangian particles (black dots), to¬ 
gether with the density field obtained from the lattice method (colored back¬ 
ground, from white: low densities to red: high densities) for St ^ 0.1 at the 
same instant of time as Fig.Q^ The fluid flow was integrated using a resolution 
of 1024^ and lattice simulations were performed with 1024^ spatial gridpoints 
associated to 17^ discrete velocities. 

To get more quantitative informations on the relevance of 
the method, we have performed a set of simulations with a 
512^ resolution and in which both the number of discrete ve¬ 
locities and the maximum velocity Fmax are varied. Fig¬ 
ure [13] shows measurements of the second-order moment of 
the coarse-grained density pr obtained by integrating the phase- 
space density /(jc, v, t) with respect to velocities and over space 
in boxes of length r. This is the two-dimensional version of 
Fig. and the statistics of (p^) have a very similar behavior as 
in the one-dimensional case. Here St = 0.5, V^ax = 3.9wnns5 
and Ny is varied from 9 to 21. One clearly observes that the 
statistics obtained from the lattice method converges to that ob¬ 
tained from Lagrangian simulations. 

The interplay between the choices of and of Vmax requires 
some further comments. On the one hand the method converges 
when both Av = 2 Vmax/A^v ^ 0 and Vmax ^ oo. On the other 
hand, the computational cost is oc One can thus wonder if 
for a flxed cost there is an optimal choice of y^ax that minimizes 



Figure 13: Second-order moment of the coarse-grained density pr as a function 
of r for St = 0.5 and TC = 1.6;: 10“Mn the direct cacade, both from Lagrangian 
measurement (black line) and the lattice method with different N^, as labeled. 

the error obtained with the lattice method. Focusing again on 
the second-order statistics of the particle mass distribution, we 
have measured the average with respect to r of the error made 
on the density moment (pj) deflned as 

S(Nv, Vmax) = <(P®)2> - <(pL)2>/<(pL)2), (14) 

where ^ 

f(r) = ^ |/(r)| rdr. (15) 



/ If rms 



Figure 14: Distance-averaged error 6 of the second-moment of the mass density 
as defined in fH) as a function of the maximal velocity Vmax (top) and of the 
velocity grid spacing Av = 2 Vmax/(A^v - 1) (bottom) for various values of the 
velocity resolution and for St ^ 0.5. 

Figure (top) represents this quantity as a function of Vmax 
for different values of the cost N^. One clearly observes that 
there is indeed for a flxed Ny a specific choice of V^ax where 
the error is minimal. The optimal value of the maximal velocity 
increases with Ny. On the right of the minimum, the error is in 
principle dominated by a Av too large. This is confirmed by the 
collapse of the various curves on the right of their minima that 
can be seen in the bottom of Fig. [T^ where 6 is represented as a 
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function of Av. In the left of the minimum, the error should be 
dominated by a too small value of Vmax- One can indeed guess 
an asymptotic collapse for Vniax ^ ^rms on the upper panel of 
Fig. or equivalently, the fact that the curves separate from 
each other at small values of Av in the lower panel. 

The value of the error at the optimal V^ax decreases from 
A^v = 9 to A^v = 13 but then seems to saturate (or to decrease 
only very slowly) at higher values of Ny. One cannot exclude 
that this behavior corresponds to a logarithmic convergence of 
S when Ny ^ oo. This slow dependence is also visible in the 
bottom panel where the collapse of the various curves seems 
to extend weakly on the left-hand side of the minima for Ny = 
13, 17 and 21. Accordingly, a small difference in Ny is not 
enough to decrease significantly the error. In the specific case 
considered (for St ^ 0.5 in the direct cascade), the resulting 
optimal choice seems to be Ay = 13 with V^ax = 2.25 Wnns. 
which leads to a relative error 10“^. 


4.2.2. Inverse energy cascade 

To complete this study we have also tested the proposed lat- 
ice method in a two-dimensional turbulent flow in the inverse 
kinetic energy cascade regime. The stochastic forcing is now 
acting at small scales (400 < |A:| < 405) and we made use of 
hyper-viscosity (here eighth power of the Laplacian) in order to 
truncate the direct enstrophy cascade. The kinetic energy accu¬ 
mulated at large scales is again removed using a linear friction 
in the Navier-Stokes equation ( p^ . The particle Stokes num¬ 
ber is now defined as St = t^Itl using the large-eddy turnover 
time tl = T/t^nns since small-scale statistics are dominated by 
forcing and are thus irrelevant. The flow is integrated with a 
resolution of 2048^ gridpoints while the lattice-particle method 
is applied for 0.1 on a coarser grid with 512^ points. 



Figure 15: Snapshot of the position of Lagrangian particles (black dots), to¬ 
gether with the density field obtained from the lattice method (colored back¬ 
ground, from white: low densities to red: high densities) for Xr ^ 0.1 in the in¬ 
verse energy cascade. The fluid flow was integrated using a resolution of 2048^ 
and lattice simulations were performed with 512^ spatial gridpoints associated 
to 9^ discrete velocities. 


Figure [T^ shows that the lattice-particle method is able to 
reproduce the main qualitative features of the particle spatial 
distribution at scales within the inertial range of the inverse en¬ 
ergy cascade. This is confirmed in Fig, which represents the 
relative error 6 defined in Eq. ( p^ of the second-order moment 
(p^) of the density pr coarse-grained over a scale r. The La¬ 
grangian integration was performed with 2x10^ particles (with 
no physical diffusion) and the lattice method on a 512^ spatial 
grid with N'l = 9^ discrete values of the particle velocity. One 
clearly observes that the error decreases at the largest scales of 
the flow. 



Figure 16: Relative error & of the second-order moment of the coarse-grained 
density pr as a function of r for St = 0.5. The lattice-particle method was here 
used with 512^ position gridpoints and 9^ velocity gridpoints. 


5. Conclusions 


We have presented a new Eulerian numerical method to 
model the dynamics of inertial particles suspended in unsteady 
flow. This lattice-particle method is based on the discretization 
in the position-velocity phase space of the evolution equation 
for the particle distribution. The spatial grid is chosen such 
that particles with a given discrete velocity hop by an integer 
number of gridpoints during one time step, an idea close to that 
used in lattice-Boltzmann schemes. We have shown that the 
model reproduces the correct dynamical and statistical proper¬ 
ties of the particles, even with a reasonably small amount of ve¬ 
locity gridpoints. Some deviations from Lagrangian measure¬ 
ments are nevertheless observed at small scales. We obtained 
evidence that they are due to numerical diffusivity acting in the 
space of velocities. The proposed method is anyway intended to 
describe large scales where such deviations disappear. It might 
hence be a suitable candidate for developing large-eddy mod¬ 
els for particle dynamics. Indeed, as equation © is linear in /, 
some techniques of subgrid modeling used in scalar turbulent 
transport (see, for example, |Girimaji and Zhou| |1996| ) could be 
adapted. 

Our approach consists in always imposing the same mesh 
for particle velocities, independently of the spatial position and 
of the local value of the fluid velocity. This is well-adapted 
for particles with a large Stokes number. Their velocity ex¬ 
periences small fluctuations and is generally poorly connected 
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to that of the fluid. In addition, the method is accurate at the 
largest scales and can hence catch the structures appearing in 
the spatial and velocity distributions of large-Stokes-number 
particles. Such considerations indicate that the proposed lattice- 
particle method is suitable for modeling particles with a suffi¬ 
ciently strong inertia. Conversely, particles with a weak inertia 
develop fine-scale structures in their distribution. They result 
from tiny departures of their velocity from that of the fluid. Our 
method, applied with a fixed velocity resolution, might not be 
able to catch such deviations. However, a more suitable idea 
for this case is to use a variation of our approach where, instead 
of a full resolution of the particle velocity, one considers its 
difference with that of the fluid. This would of course require 
changing scheme for integrating advection. 
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